as = [];
for alpha=0.01:0.01:0.99
    h=0.01;N=100;z=sqrt(3)/pi*(log(alpha)-log(1-alpha));
    y(1)=h;y(2)=h+h*(h*(1+z)+1);
    for i=2:N
        y(i+1)=y(i-1)+h/3*(y(i)*(1+z)+1)+1/2*h^2*(y(i)*(1+z)+1)*(1+z)-1/6*h^2*(y(i-1)*(1+z)+1)*(1+z)+4/3*h*(y(i)*(1+z)+1)+h/3*(y(i-1)*(1+z)+1);
    end
    as = [as;y(i)];
end
plot(as,0.01:0.01:0.99,'g-.')
